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We report diffusion quantum Monte Carlo calculations of three-dimensional Wigner crystals in the 
density range r 3 — 100 to 150. We have tested different types of orbital for use in the approximate 
wave functions but none improve upon the simple Gaussian form. The Gaussian exponents are 
optimized by directly minimizing the diffusion quantum Monte Carlo energy. We have carefully 
investigated and sought to minimize the potential biases in our Monte Carlo results. We conclude 
that the uniform electron gas undergoes a transition from a ferromagnetic fluid to a body-centered 
cubic Wigner crystal at r 3 = 106 ±1. The diffusion quantum Monte Carlo results are compared with 
those from Hartree-Fock and Hartree theory in order to understand the role played by exchange and 
correlation in Wigner crystals. We also study "floating" Wigner crystals and give results for their 
pair correlation functions. 

j3 1 PACS numbers: 71.10.Ca, 71.15.Nc, 73.20.Qt 
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I. INTRODUCTION 



D | Electrons in a uniform potential are expected to crystallize at low densities to minimize their interaction energy* 1 * 2 - 
;_i ■ In this paper we investigate Wigner crystals in three dimensions using quantum Monte Carlo (QMC) methods within 
the variational (VMC) and diffusion (DMC) approaches. The DMC method is currently the most accurate available 
for calculating the zero-temperature ground state energy of extended quantum mechanical systems. 

There has been some debate about the density at which the transition from a Fermi fluid to a Wigner crystal 
should occur in three dimensions. In their pioneering DMC study of the phases of the electron gas, Ceperley and 
Alder 3 obtained a transition density of r s = 100 ± 20, 4 but a more recent DMC study gave r s — 65 ± 10£ Moreover, 
highly accurate DMC energies for the low-density fluid have recently become available^, which may further modify 
predictions of the transition density. The primary goal of this work is to provide highly accurate DMC energies for 
£j \ three-dimensional Wigner crystals and to use them in conjunction with the fluid data of Zong et alJ± to predict an 
accurate value for the transition density. 

To achieve sufficient accuracy we have carefully studied the possible sources of error in our calculations, including 
finite size effects, timestep errors and population control errors. We have used three procedures for optimizing the 
' trial wave functions: minimization of the variance of the energy within VMC-' ? ; minimization of the VMC energy; 
, and minimization of the DMC energy. 

We compare our results with recently-published fully-self-consistent unrestricted Hartree-Fock (HF) calculations 9 , 
and with the results of a simple version of Hartree theory, which allows us to understand the effects of exchange and 
— i . correlation in Wigner crystals. 

Finally, we discuss "floating" Wigner crystals in which the homogeneous and isotropic nature of the ground state is 
restored, relating their properties to those of a "fixed" crystal. DMC results for the pair correlation functions (PCFs) 
' of floating Wigner crystals are compared with those of the fluid phases. 



II. QUANTUM MONTE CARLO METHODS 

We have performed both VMC and DMC calculations using the CASINO packaged In VMC, expectation values are 
calculated using an approximate many-body wave function, the integrals being performed by a Monte Carlo method. 
The approximate wave function normally contains a number of variable parameters, whose values are obtained by an 
optimization procedure. 

In the DMC method^ 1 - the imaginary time Schrodinger equation is used to evolve an ensemble of electronic 
configurations towards the ground state. The fermionic symmetry is maintained by the fixed-node approximation 12 , 
in which the nodal surface of the wave function is constrained to equal that of an approximate wave function. We 
will refer to the approximate wave functions used in VMC and DMC as trial wave functions. The fixed-node DMC 
energy provides a variational upper bound on the ground state energy with an error that is second order in the error 
in the nodal surface i 13 ! 14 

The trial wave function introduces importance sampling and controls both the statistical efficiency and the final 
accuracy that can be obtained. In DMC the final accuracy depends on the nodal surface of the trial wave function 
via the fixed-node approximation, while in VMC the final accuracy depends on the entire trial wave function, so 
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that VMC energies are more sensitive to the quality of the trial wave function than DMC energies. Apart from the 
fixed-node approximation, DMC results may be subject to bias from the use of the short-time approximation (finite 
timestep errors), population control errors, and effects from the finite size of the simulation cell. We have made 
strenuous attempts to reduce these errors: see Sec. IIVI The statistical errors in our QMC data are estimated using 
the blocking method^ to eliminate the effects of serial correlation. 

III. TRIAL WAVE FUNCTIONS 

A. The Slater-Jastrow form 

We have used trial wave functions of the standard Slater-Jastrow form, 

* (R T , RJ = e'C^'^Dt (R T ) D l (R x ) , (1) 

where and are Slater determinants of up and down spin single-particle orbitals and Rj and Rj denote the 
coordinates of the up and down spin electrons, respectively, and e J is the Jastrow correlation factor. 

B. Jastrow factors 

We have used Jastrow factors of the form 

i 3 

where 

Mrv) -±(l-„(-%))~,(-$), (3) 

with Fij = y2A if electrons i and j have the same spin and F^ = \fA if the electrons have opposite spins. This 
term satisfies the electron-electron cusp conditions . 11 i 16 The constant Lq is set to 0.3 of the Wigner-Seitz radius of 
the simulation cell and A is a free parameter. The uq term is set to zero for greater than the Wigner-Seitz radius, 
resulting in a small discontinuity in the Jastrow factor of less than 2 x 10~ 5 in magnitude. To investigate possible 
bias from this discontinuity we compared the values of the two standard estimators of the kinetic energy involving 
the gradient and Laplacian of the trial wave function from a very long VMC run at r s = f 00. The estimators agreed 
to within the statistical error, which was smaller than in our final DMC runs. 
The second term in the Jastrow factor is given by 

L 

+B'(r ij -L'f(^- + r ij y (4) 

where Tj is the Ith Chebyshev polynomial, L' is the Wigner-Seitz radius of the simulation cell and B' and the ai are 
parameters to be determined. 

C. Orbitals for the Slater determinants 

The Slater determinants for the crystalline phases were formed from localized non-orthogonal single-particle orbitals 
centered on the lattice sites of a body-centered cubic (BCC) crystal. A BCC crystal is expected to be favored in the 
low-density limit because it has the lowest Madelung energy. Throughout, we use 0(r) to denote a spatial orbital 
centered on the origin. Periodic orbitals for use in a simulation of a finite system subject to periodic boundary 
conditions are constructed for each lattice point in the simulation cell by summing over all the replicas of (f> centered 
on that lattice point. Clearly, if all the individual orbitals are periodic then their Slater determinant is too. We use a 
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Jastrow factor containing only homogeneous terms, see Sec. IIIIB1 with the differences between electron coordinates 
being evaluated under the minimum image convention. Hence the overall wave function is periodic. 

In their VMC and DMC studies of Wigner crystals, Ceperleyi^ and Ceperley and Alder— used Gaussian orbitals: 

0(r) = exp(-CV 2 ). (5) 

They determined C by variational methods, with the Jastrow factor being optimized simultaneously. 
Ortiz et aZ.~ used exponentials of two-parameter Pade functions: 

^)=exp(^) ) (6) 

and determined the values of C and k and the parameters in their Jastrow factor by minimizing the variance of the 
energy within VMC. 

We have also investigated two new types of orbital for Wigner crystals. A straightforward generalization of Eq. (O 
is to use a linear combination of Gaussian orbitals: 

N G 

0(r) =^A l exp(-C,r 2 ). (7) 

i=i 

We have also considered orbitals based on an expansion in the eigenstates of a simple harmonic oscillator. An orbital 
is constructed by multiplying the simple Gaussian function by a polynomial. For a BCC lattice with identical orbitals 
on each site the polynomial should have the full symmetry of the lattice, i.e., 

4>{y) = cxp(-CV 2 ) [l + ar 2 + f3r A + j(x 2 y 2 + x 2 z 2 + y 2 z 2 ) + 0(r e )] , (8) 

where r = (x,y, z). This orbital has considerable flexibility at small r. 

D. Optimization of the trial wave functions 

Parameters in the trial wave functions may be optimized in a variety of ways. In principle the DMC energy depends 
only on the nodal surface of the wave function, which is determined by the form of the orbital. It is therefore best to 
minimize the DMC energy directly with respect to the parameters in the orbitals, but this is a costly and laborious 
procedure which we have carried out only for the Gaussian parameter C of the simple Gaussian orbital. In principle 
the DMC energy does not depend on the Jastrow factor, so it cannot be optimized in this fashion. 

We first studied the simple Gaussian and Pade forms of Eqs. (JSJ) and ©, using energy variance minimization to 
optimize the C parameter and the parameters in Jastrow factor, and minimization of the VMC energy with respect 
to variations in k, but we found the optimal value of k to be very close to zero. The Pade orbital seems to offer little 
advantage over a simple Gaussian orbital at the densities studied (100 < r s < 150). 

We optimized the expansion coefficients and Gaussian exponents of the linear combination of Gaussian orbitals 
form of Eq. ([7]), together with a Jastrow factor, using VMC energy variance minimization. However, again, it was 
found that in practice this orbital offers no advantage over a single Gaussian function at the densities studied. 

We attempted to optimize the a parameter in the harmonic oscillator form of Eq. ([5]) for a Wigner crystal at 
r s = 100, but the resulting wave function gave a DMC energy within the error bars of the one obtained by setting 
a = 0. 

We have therefore adopted the simple Gaussian orbital of Eq. ([5]) for our main calculations. We have adopted the 
following procedure to optimize the trial wave functions. The A parameter in the Jastrow factor was optimized by 
minimization of the VMC energy, the parameters in the Si part of the Jastrow factor by energy variance minimization, 
and the C parameter in the Gaussian orbitals by minimization of the DMC energy. These minimizations were 
performed in turn until the changes in the parameters were negligible. We found that the DMC-optimized exponents 
obey C » 0.11r7 3 ^ 2 , which gives rather smaller values than those used by Ceperley2iIZ (C « 0.2r s 3 ^ 2 ) and considerably 
smaller than those predicted by HF theory or the simple Hartree model (C ~ 0.5r s 3 ^ 2 ), see Sec. IVII1 

IV. ACCURACY OF THE DMC RESULTS 

Our DMC algorithm is essentially that of Umrigar et al*& Here we explore the sources of error in our DMC 
calculations and justify our choices of the parameters for the final production runs, which are summarized in Sec. lIVDl 
Unless otherwise stated, we use simple Gaussian orbitals throughout this section. 
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A. Finite size effects 

We used periodic boundary conditions and the Ewald interaction energy to reduce the finite size effects. We tested 
the convergence of the Ewald sums and found that truncation errors were less than 10 -3 of the statistical error bars 
on the final DMC runs. 

The energy per electron at a given density depends on the number of electrons in the simulation cell. We wish 
to obtain the energy per electron in the limit that the number of electrons per simulation cell goes to infinity. Two 
approaches have been used previously when dealing with finite size effects in QMC simulations of Wigner crystals. 
Ortiz et alM used large simulation cells and found the finite size errors to be less than their statistical error bars of 
8.5 x 10~ 6 a.u. per electron for numbers of electrons in excess of 500. Ceperley 3 -^, on the other hand, used smaller 
system sizes in conjunction with an extrapolation scheme. 

Because we wish to work to very high accuracy, we use quite large simulations cells and the extrapolation formula 
derived by Ceperleyii 

b 

Eoo = En H — 575 — , (9) 
r7 2 N 

where the constants Eoo and En are the total energies per electron of the infinite system and the system with N 
electrons, and b is independent of both r s and N. 

Starting from a two-parameter Pade orbital and corresponding Jastrow factor optimized in a 64-electron simulation 
cell, we attempted to further optimize the wave function in a 216-electron unit cell. This attempt did not lead to a 
lowering of the VMC energy, suggesting that the 64-electron simulation cell is adequate for optimization purposes; 
this size of cell was used in all of our subsequent optimization runs. 

B. Population control biasing 

The use of a finite population of configurations results in a positive bias in the DMC energy which, it is argued, 
falls off as the reciprocal of the target population. 18 This turns out to be a genuine problem in the case of Wigner 
crystals where we are able to work to extremely high precision. An example of the problem of population control 
biasing is shown in Fig. [TJ where it can be seen that the bias is indeed positive and that it falls off roughly as the 
reciprocal of the target population. 
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FIG. 1: DMC energy against target population for a 64-electron crystal at r 3 — 100. The Gaussian exponent was C = 0.000135 
and the Jastrow factor contained only the uq term, with A = 438.389. The timestep was 20 a.u. 

The simplest method for avoiding population control biasing is to use a large target population. Alternatively, 
the reweighting scheme developed by Umrigar et alJ£ can be used. In our tests we found this scheme to work well, 
provided the number of reweighting factors was about the same as the number of timesteps over which average local 
energies are correlated (between 100 and 1000 timesteps of 30 a.u.). However, with this many reweighting factors 
present, the total weight at each timestep fluctuated enormously and very long simulations were required in order to 
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obtain meaningful statistics. We found it to be more efficient to use larger populations than to employ the reweighting 
scheme. For this reason, we did not use the reweighting scheme in our production runs. 

Including more parameters in the Jastrow factor can lead to a significant reduction in population control biasing, 
as illustrated in Fig. [Jl This shows that, when the reweighting scheme is not used, DMC energies obtained with a 
a poor Jastrow factor and a small target population (solid line) are too high, but when a larger population is used 
in conjunction with the same Jastrow factor (dashed line) the energies are similar to those obtained with a good 
Jastrow factor and a small population (dotted line). Improving the overall quality of the trial wave function reduces 
the population control bias because it reduces the fluctuations in the population. The results shown in Fig. [1] in which 
the Jastrow factor consisted only of the uq term therefore represent a worst-case scenario. 
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FIG. 2: DMC energy against Gaussian exponent for a 64-electron crystal at r s — 125, using a timestep of 30 a.u. and A — 
597.901. Solid line: DMC energies with a Jastrow factor consisting of only the wo term and a target population of 100 
configurations; dashed line: DMC energies for the same Jastrow factor, but with a target population of 800 configurations; 
dotted line: DMC energies obtained with a fully optimized Jastrow factor (with both uo and Si terms) and a target population 
of 100 configurations. 



C. Timestep biasing 

The variation of DMC energy with timestep is shown in Fig. [3] for three different Jastrow factors. It can be seen 
that the bias is always positive and that it grows linearly with timestep; therefore we can largely eliminate the bias by 
carrying out simulations at a number of different timesteps and performing a linear extrapolation to zero timestep. 

The differences between the DMC energies in Fig. [3] are due to population control and timestep bias. The solid 
and dashed curves were obtained using the same target population of 100 configurations, but with different Jastrow 
factors. The population control bias at fixed timestep is clearly smaller for the Jastrow factor with A — 438.389 (solid 
line) than for A — 563.157 (dashed line). The dotted line was obtained with a target population of 800 configurations, 
which essentially removes the population control bias. Because the solid and dotted curves are approximately parallel 
we deduce that the population and timestep errors are approximately independent of one another. Furthermore, it 
is clear that altering the Jastrow factor has a considerably greater effect on the population control bias than on the 
timestep bias. 

In Fig. [4] we show that timestep bias remains a problem even with a well-optimized Jastrow factor and a large 
target population which essentially removes the population bias. However, this figure again shows that a linear fit is 
appropriate when extrapolating to zero timestep. 

D. Parameters for the production runs 

The final production runs were characterized as follows. 

1. The Gaussian exponents in the orbitals were optimized by minimizing the DMC energy, the A parameters by 
minimizing the VMC energy, and the other parameters in the Jastrow factors by minimizing the variance of the 
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FIG. 3: DMC energy of a 64-electron crystal at r s = 100 against timestep for different values of the A-parameter in the Jastrow 
factor of Eq. © and different target population sizes. The Si term was not present in the Jastrow factor. The Gaussian 
exponent was C = 0.00011 in all cases. Solid line: A — 438.389, target population 100 configurations; dotted line: A — 438.389, 
target population 800 configurations; dashed line: A — 563.157, target population 100 configurations. 
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FIG. 4: DMC energy of a 64-electron crystal at r s — 110 against timestep. The Jastrow factor contained optimized uo and Si 
terms and the target population was 640 configurations. The straight line is a fit to the DMC data. 



energy, as described in Sec. IIII Dl The Si terms in the Jastrow factors contained between 4 and 6 parameters 
per spin. 

2. A target population of 640 configurations was used. This, together with the optimized Jastrow factor, should 
ensure that population control errors are negligible. 

3. At each density, DMC calculations were performed using between four and six different timesteps and the energy 
was extrapolated linearly to zero timestep. 

4. A variety of system sizes were used (see Table HT]), and the energies were extrapolated to infinite system size 
using Eq. ©. 
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V. RESULTS AND DISCUSSION 
A. DMC energies obtained using DMC-optimized Gaussian orbitals 

The values of the exponents obtained by optimizing the DMC energy are shown in Table HI along with the final 
DMC energies. These were found by using the DMC results shown in Table |TT] in conjunction with the finite size 
extrapolation formula of Eq. ([9]). Using the results for r s = 100 and r s — 125, we obtained a good fit, giving 
6=1.26(3). 



r 3 


C (DMC) 


DMC energy 


100 


0.00011 


-0.0076765(4) 


110 


0.0001 


-0.0070312(5) 


125 


0.00009 


-0.0062458(4) 


150 


0.000063 


-0.0052690(3) 



TABLE I: Orbital exponents optimized by minimizing the DMC energy, and the final DMC energies per electron (extrapolated 
to zero timestep and infinite system size) for the different r s . All entries are in a.u. 



r s 


System size 


DMC energy 


100 


64 


-0.0076961(2) 


100 


216 


-0.0076823(3) 


100 


512 


-0.0076788(8) 


110 


64 


-0.0070483(2) 


125 


64 


-0.0062599(2) 


125 


216 


-0.0062495(1) 


150 


64 


-0.0052797(1) 



TABLE II: DMC energies in a.u. per electron (extrapolated to zero timestep) at different densities and system sizes, which are 
characterized by the number of electrons in the simulation cell. 



B. Electronic charge densities obtained using the different orbitals 

In Fig. [5] we plot the electronic charge densities for a r s = 100 crystal, calculated using HF theory orbitals (see 
Sec. IVTip . from the DMC-optimized orbitals but without a Jastrow factor, and within DMC using an optimized 
Jastrow factor. 

The HF theory orbitals are very localized whereas the orbitals optimized within DMC are much more diffuse. 
However, the inclusion of a Jastrow factor results in the charge density from the DMC-optimized orbitals becoming 
more localized on lattice sites, although not to the same extent as the HF orbitals. The VMC charge density obtained 
with the optimized Slater- Jastrow wave function is very close to the DMC density shown in Fig.O The peak difference 
between the DMC and VMC charge densities is 2.6 arb. units, and therefore an extrapolated estimation of the QMC 
charge density, p{r) ss 2prjMC ~ Pvmc> would be very similar to the DMC density. 

The Jastrow factor in a Wigner crystal serves to further localize the electrons on their lattice sites. Therefore, 
whereas HF theory gives very localized Gaussian orbitals, in an optimized Slater-Jastrow wave function we find the 
orbitals to be much more diffuse, with the localization being caused by correlation effects from the Jastrow factor. 
Correlation effects allow electrons to invade each other's space to some extent without incurring a high potential 
energy penalty; hence the overall charge density is less localized than in HF theory, and the kinetic energy is lower as 
a result. 
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FIG. 5: Electronic charge density along a (111) direction passing through the lattice sites of a BCC Wigner crystal at r 3 = 100. 
Solid line: HF charge density; dotted line: charge density from the DMC-optimized orbitals without a Jastrow factor; dashed 
line: DMC charge density. 



VI. LOCATING THE FLUID-TO-CRYSTAL TRANSITION DENSITY 



In order to locate the density at which the transition from the Fermi fluid to the crystal occurs, we fit the DMC 
energy data for these phases to interpolating functions. The correlation energy is defined as 

E c (r s ,0 = E(r a ,0-E HF (r s ,0, (10) 

where E is the total energy, Ehf is the HF ground state energy, £ is the polarization. Ceperleyil proposed a fitting 
form for the correlation energy of a Fermi fluid as a function of r s , 

E c( r s) = - ~ C Ti ~ C , (11) 

where 7^, /?£ and /?£ are fitting parameters. We make use of the highly accurate DMC energies of Zong et al& for the 
ferromagnetic fluid phase at low densities, which used trial wave functions including "backflow" effects. We found an 
excellent fit to Eq. {nj giving 71 = -0.09399, 0{ = 1.5268 and 0\ = 0.28882. We also tried fitting the fluid data to 
the form proposed by Perdew and Zunger 19 , which is based upon Eq. |[TT|) . but includes an assumed dependence on 
polarization, so that the partially polarized and unpolarized data of Zong et al. could be used. The \ 2 value of this 
fit was not so good, however. 

At low densities the total energy of a crystal phase can be expanded as 

^) = -+^+4+<^ 5/2 )> (12) 

r s r */ z r% 

where the {/} are constants^ The first of these is taken to be /o = —0.89593 in order to give the Madelung energy in 
the low density limit. We found that our DMC data fitted Eq. (O very well, giving f x = 1.3379 and f 2 = -0.55270. 
These values are in reasonable agreement with the parameters found using a completely different method by Carr— 
and Carr, Coldwell-Horsfall and FeinSI, who have calculated the zero-point lattice-vibrational energy of a Wigner 
crystal in order to give an analytical result of f% = 1.325. Furthermore, they use perturbation theory to obtain the 
approximate result f 2 = —0.365. This phonon model is in good agreement with our Wigner crystal energies at large 
r s , but it gives energies which are too high at smaller r s . 

The energies of the ferromagnetic fluid and BCC crystalline phases at low densities calculated by different authors 
are shown in Fig. [6] We found the transition from the fluid to crystalline phases to occur at r s = 106 ± 1, in agreement 
with the original result of Ceperley and Alder^ Note, however, that the transition density predicted using the fluid 
data of Zong et al. 6 in conjunction with the crystal data of Ceperley and Alder would be somewhat lower, at about 
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r s = 127. Our Wigner crystal energies are slightly lower than those of Ceperley and Alder, even though they studied 
a Bosc crystal, which must have a lower energy than the corresponding fcrmion one. We believe this difference must 
be due to some small bias in the results of Ceperley and Alder.— Fitting Eq. (|12jl to the crystal data of Ceperley 
and Alder, we find that fx = 1.4309 and = —1.1058. The discrepancy with the analytical result of Carr for f\ is 
consistent with the presence of a small, positive, systematic bias in the crystal results of Ceperley and Alder. 

Our transition density is considerably lower than the value of r s = 65 ± 10 obtained by Ortiz et aL— >22. The 
statistical error bars on their data are much larger than on the fluid data of Zong et al& or on our Wigner crystal 
data, which hampers detailed comparisons. However, it appears that the main reason for the discrepancy is that 
Ortiz et al. place the ferromagnetic fluid higher in energy than Zong et al. Some bias in the results of Ortiz et al. is 
expected because they used plane wave nodes while Zong et al. used optimized backflow nodes, but this is not enough 
to explain the difference between the results. 
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FIG. 6: Energies of the ferromagnetic fluid and BCC crystalline phases at low density. The Madelung energy of the BCC lattice 
has been subtracted off and the resulting energy multiplied by r^ 2 to highlight the differences between phases. The circles 
show the DMC data of Zong et al& for the ferromagnetic Fermi fluid; the dashed line is a fit to this data. The diamonds show 
our DMC results for the BCC crystal; the solid line is a fit to this data. The left-pointing triangles show the Ceperley-Alder- 
results for the ferromagnetic fluid; the dotted line is a fit to this data. The up-pointing triangles show the Ceperley- Alder results 
for the BCC crystal; the dashed-dashed-dotted line is a fit to this data. Finally, the dashed-dotted line shows the prediction of 
the phonon model of Carr et al?- Where error bars on the DMC data cannot be seen, it is because they are smaller than the 
symbols showing the data points. 

For the crystal phase, we may estimate the fixed node error resulting from the use of the HF orbitals by taking 
the difference between DMC energies calculated using the HF and DMC-optimized orbitals. At r s — 100 we find this 
difference to be 8.9(1) x 10 -5 a.u. per electron. Zong et al~2- have calculated the fixed node errors resulting from 
the use of plane-wave orbitals for the fluid phases as 1.87(5) x 10 -5 a.u. per electron for the unpolarized fluid and 
0.84(7) x 10 -5 a.u. per electron for the fully polarized fluid. Therefore, although the correlation energy of crystals is 
smaller than fluids at the same density, the fixed node errors resulting from the use of HF orbitals are considerably 
larger in crystals than in fluids. 

VII. HF AND OTHER SIMPLE THEORIES 

HF theory is described as "restricted" when the spin-orbitals are products of space and spin parts which are 
occupied in pairs with identical space parts, and "unrestricted" when the space parts are different or when they are 
not occupied in pairs. Within the quantum chemistry community the standard definition of the correlation energy 
is the difference between the exact and restricted HF energies, but in a Wigner crystal the electrons are localized in 
space individually and a description within restricted HF theory is not possible. We therefore define the correlation 
energy as the difference between the exact and unrestricted HF energies. 

A recent self-consistent unrestricted HF study of electrons in a uniform potential gave stable Wigner crystal solutions 
for r s > 4.5 in three dimensions.— Here we develop simplified versions of the HF model which almost exactly reproduce 
the results of the fully self-consistent HF studies^ at low densities and compare the results with our DMC ones. 
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In the low-density regime, the overlap between orbitals centered on neighboring lattice sites is small and therefore 
we expect the Hartree and HF energies to be similar. Let us take the wave function of the crystalline phase to be a 
Hartree product of normalized Gaussian orbitals, 

2cX 3/4 



<H r ) = [—) exp(-CV), (13) 

centered on lattice sites. The Gaussian exponent C is to be determined variationally. 

The resulting kinetic energy per electron is easily evaluated as T — 3C/2. The spatial charge density is simply 
the superposition of the Gaussian charge densities due to each orbital. The electrostatic energy Eh of this charge 
distribution may be readily evaluated, but we must subtract off the self-energies Eq of each Gaussian. The total 
energy per electron is therefore given by 

E = T + Eh — Eq 

_ 3C 2tt v cxp(-G 2 /4C) [C 

2 n ^ g 2 V 7T ' [ ' 

G^O 

where the G are the reciprocal lattice vectors and f2 = 4irr^/3 is the volume of the primitive unit cell. 

Differentiating this energy with respect to the Gaussian exponent and approximating the sum by an integral, we 
find that 

(15) 

dC 2 2SIC 2 ' V ; 

The approximation is valid for large r s , where the density of reciprocal lattice vectors is large. Hence we find 
the optimal value of C to be C = l/2r^ 2 , which is precisely the result obtained by Wigner— using a spherical 
approximation. 

Integrating Eq. (|15p with respect to C we obtain the energy E = 3C/2 + ir/2£lC + f(r s ), where the function f(r s ) 
is the "constant" of integration. Inserting the optimal value of C and making use of the fact that, in the limit of low 
densities, E must tend to the Madelung energy of the crystal lattice, we find 

E =Aj2+-' ( 16 ) 

where M is the Madelung constant of the lattice. 

The simple Hartree model of Eq. (| 14[) and the further simplifying approximation of Eq. (|16|) give energies which 
are very close to the full HF results for r s > 50; for example, the energy of Eq. (TTJJ agrees with that of the full 
HF result to within 0.006%, whereas the energies at r s — 100 agree to within 0.001%. The agreement between 
the fully self-consistent HF data and the simpler approximations is just as good for a face centered lattice. These 
results demonstrate that exchange energies between orbitals on different sites are extremely small. The exchange 
interactions are only significant between nearest neighbor Gaussian orbitals and the exchange energy per electron 
is given by Ex = —(N n /2)\/C/irex.p(—CR 2 ), where N n is the number of nearest neighbors and R is the nearest 
neighbor distance. This expression gives extremely small energies for the low densities studied here. 

We have also calculated DMC energies using the orbitals obtained from HF theory^, and Jastrow factors optimized 
using energy variance minimization. The results are summarized in Table [THl The extrapolation of the DMC results 
to infinite system size was carried out using Ceperley's extrapolation scheme with his value of b = I.5. 1 - Comparing 
the data in Table [TTT] with that in Tables [J or HQ we see that the DMC energy obtained at r s — 100 with the HF 
orbitals is significantly higher than that obtained with DMC optimized orbitals. Some population control bias is 
present in the results in Table HTT1 but we estimate this to be about 2% of the difference between the r s = 100 DMC 
energies obtained with the HF and DMC-optimizcd orbitals. 



r s 


C (HF) 


HF energy 


DMC 


energy 








(512 electrons) 


(Infinite system) 


50 


0.00141 


-0.0136768 


-0.014060(4) 


-0.014052(4) 


100 


0.0005 


-0.0074593 


-0.007589(1) 


-0.007586(1) 



TABLE III: HF results and DMC results obtained using HF theory orbitals. All entries are in a.u. 
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The strength of correlations in a system may be measured by the ratio of the correlation energy to the total energy, 
E c /E. The DMC results of Zong et al£ indicate that in the fluid phases E c /E tends to a positive constant as r s — ► oo, 
while for the Wigner crystal our results show that E c /E tends to zero as r s — > oo. In this sense one may think of 
Wigner crystals as being weakly correlated systems at low densities. 

VIII. MAGNETIC BEHAVIOR OF THE CRYSTALLINE PHASES 

The tiny energy differences between ferromagnetic and antiferromagnetic crystals proved to be too small to resolve 
in our QMC calculations. It might be possible to resolve them using a correlated sampling approach within VMC. Such 
an approach should provide an accurate value for the energy difference between two systems, 1 and 2, if \^i\ 2 ~ l^l 2 
throughout the configuration space, which should hold for ferromagnetic and antiferromagnetic crystalline phases at 
sufficiently low densities. HF theory predicts ferromagnetic behavior in the low density limit but according to the 
theory of Thouless- 2 ^ such a system should be antiferromagnetic. In their path integral QMC calculations, Candido, 
Bernu, and Ceperley 2 ^ have indeed found antiferromagnetic behavior for BCC Wigner crystals at low densities, 
although the energy differences are much smaller than our statistical error bars. 

IX. FLOATING WIGNER CRYSTALS 

The Hamiltonian of the uniform electron gas is invariant under the simultaneous translation or rotation of the 
electron positions. However, our trial Wigner crystal wave functions break these symmetries and, for example, the 
resulting charge densities are inhomogeneous, see Fig. [5] These wave functions represent Wigner crystals which are 
"pinned" by some small external influence. Pinning of Wigner crystals may arise from the presence of impurities or 
boundaries, and therefore the broken symmetry solutions are physically relevant, but Wigner crystals may also be 
mobile, in which case it is better to describe them as floating crystals. 

One way of restoring the homogeneous and isotropic nature of the trial function is to consider a linear combination 
of displaced and rotated copies of the fixed wave function \F This gives a "floating" wave function, 

1 fF = J^{R{^){{r l -A}))dndA, (17) 

where R((l) represents a rotation of all of the electron positions and the integrals are over all possible solid angles 
17 and displacements A. gives rise to a homogeneous and isotropic charge density n(r) = N/Sl, where N is the 
number of electrons and fl is the volume of the system. 

Wave functions for floating Wigner crystals in extended systems have been discussed briefly by, for example, 
Bishop and Liihrmann^, and by Mikhailov and Ziegleri^l Rather more attention has been devoted to restoring the 
(rotational) symmetry in two-dimensional models of quantum dotsi 27 ' 28 In finite systems the energy gain per electron 
from restoring the symmetry can be substantial but for an infinite system it turns out to be negligible. 

Using a trial function where the single particle orbitals in the Slater determinant are <j> — e~~ Cr , we have obtained 
analytical results at the variational level for the case when the translational symmetry is restored. We found that 
the difference in total energy is equal to the kinetic energy of the center of mass of the fixed crystal (3C/2), making 
the energy difference per electron negligible (this result also holds when a translationally invariant Jastrow factor is 
included) . We also found that the expectation value of any operator that only depends on relative electron coordinates 
is the same for >F and $p at the variational level. We expect the same conclusions to be true in DMC as well; in 
particular, we found that under open boundary conditions the nodal surface of the fixed and the translationally 
averaged trial functions are identical. 

Despite the similarities in their energies, there is an important qualitative difference between the fixed and floating 
Wigner crystal wave functions. The fixed wave function can be written as a sum of disconnected partial wave functions, 
in the form $ = X^mV'm, where the overlap between the -0m tend exponentially to zero as N — ► oo. From this it 
follows that it represents an insulator i 29 i 30 On the other hand, the same is not true for the floating wave function, 
resulting in conducting behavior. 
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X. PAIR CORRELATION FUNCTIONS 
A. Definitions 

The spin-dependent PCF is defined as 



g*A r y) = - — — n — tt, -• ( 18 ) 

It would be very costly to evaluate the full six-dimensional function <7 CT)CT ' (r, r') for a Wigner crystal within QMC. 
However, for a homogeneous and isotropic system g depends only on the separation between electrons, r = |r — r'|, 
so that 



^r) = ^ 7 ^(^Sir-\r i -r j \)Y (19) 



where we have used n a (r) — N a /Vl. This one-dimensional function is much less costly to evaluate accurately than 
Eq. (|18p . hence we can obtain the PCF for a floating Wigner crystal. Since the expectation value of any operator that 
depends only on the difference between electron coordinates is the same for fixed and floating crystals, see Sec. IIX| 
we can evaluate g a ,a' (f) using our fixed trial function. 

Furthermore, the expectation value of the potential energy of the system can be written in terms of the spin- 
independent PCF, g(r) = jJ2a,a> 9<y,a'{r): 



J^v^ri-rjDj = J n{v)n{v')g{v,v')v{\v-v'\)dvdv' 



n 2 r 

= — / 4nr 2 g(r)v(r)dr, (20) 

which holds even for an inhomogeneous system such as the fixed Wigner crystal. Studying the PCF of a floating 
crystal therefore provides insight into the physics of the fixed crystal as well. 



B. Discussion and Results for the PCFs 



We evaluated Eq. (TT9^) within VMC and DMC by accumulating g a a ,(r) in radial bins. Our best estimates of g 
were obtained using the extrapolated estimator g[r) ~ 2^dmc ~ 5vmc- We calculated g a .cr l (j') at r s = 110 for a 
ferromagnetic and an antiferromagnetic Wigner crystal using our optimized trial wave functions. For comparison, 
we have also calculated g a , a '{r) for the unpolarizcd fluid phase at r s = 110, using a trial wave function consisting of 
determinants of plane waves and a Jastrow factor optimized using energy variance minimization. 

All of the biasing effects that apply to the DMC energy may also affect the PCFs. We found that this was indeed 
the case, as we have obtained results that showed very small but statistically significant differences when, for example, 
different timesteps were used or when different Jastrow factors were used in conjunction with small population sizes. 
Unlike the energy, it was not possible to use an extrapolation scheme to remove the timestep bias, as it showed no 
clear pattern. We found that finite size effects in the PCFs were very small, however, as results obtained for 64 
electrons were not significantly different from those obtained with 512 electrons. 

A further source of bias, which does not apply to the energy, arises from the use of the extrapolated estimator. 
Fig. [7] shows that the extrapolation can make a significant difference, and to check the reliability of this method we 
have evaluated the PCF using different quality Jastrow factors. As might be expected, the VMC and DMC results 
varied significantly, but the final extrapolated results were almost independent of the Jastrow factor, the differences 
being only slightly larger than the statistical error. This source of bias is therefore small and on a comparable level 
to the others. 

The final PCF calculations were carried out in a 512 (518) electron system for the Wigner crystal (fluid), using 
a timestep of 30 a.u. and a target population of 960. Fig. [H] shows the extrapolated spin-independent PCF for the 
antiferromagnetic crystal and unpolarized fluid phases at r s = 110. For comparison, the HF result is also shown, 
together with the result obtained from the DMC-optimized orbitals, but without the Jastrow factor. The spin- 
independent PCF for the ferromagnetic crystal is not included as it was found to be almost identical to that of the 
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FIG. 7: Spin-independent PCF in the antiferromagnetic floating crystal at r s — 110, illustrating the extrapolation procedure. 
Solid line: extrapolated estimator; dashed line: DMC result; dotted line: VMC result. The statistical errors are less than 0.003. 



antiferromagnetic crystal. The HF orbitals are very strongly localized and give the most rapidly varying PCF, whereas 
the PCF from the more diffuse DMC-optimized orbitals is much smoother. 

It is interesting that the extrapolated PCF of the floating crystal shows strong oscillations at distances much larger 
than r s . This can be understood in terms of Eq. (|20p . The potential energies of the fixed and floating crystals are the 
same, but their charge densities and PCFs are quite different. For the floating crystal n(r) is constant while g(r, r') 
oscillates strongly, whereas for a fixed crystal the charge density n(r) oscillates and g(r, r') is expected to be much 
smoother. 

Another interesting point is the large difference between the extrapolated parallel/antiparallel-spin PCFs of the 
crystal and the fluid, see Fig. [9l For the crystal, the PCF strongly reflects the underlying crystal structure composed 
of alternating spins, whereas for the fluid phase the two are almost identical. The energies of the two systems are, 
however, almost the same. 

XI. CONCLUSIONS 

We have carried out a careful DMC study of the BCC Wigner crystal in the density range 100 < r s < 150. We 
have experimented with several types of orbital in the trial wave functions but have been unable to improve upon the 
Gaussian form used in previous work. We have, however, optimized the Gaussian exponent by directly minimizing 
the DMC energy, which reduces the fixed node errors. We have also taken care to eliminate other biases in our DMC 
simulations, particularly those from timestep errors, population control bias, and finite size effects. We estimate that 
the uniform electron gas undergoes a transition from a ferromagnetic fluid to a BCC crystal at r s = 106 ± 1. 

We have used Slater- Jastrow- type trial wave functions for our studies. Multiplying the Slater determinant by a 
pairwise repulsive Jastrow factor makes the charge density more inhomogeneous because the electrons in Wigner 
crystals are localized in space individually. This behavior contrasts with that found in many other systems where a 
pairwise repulsive Jastrow factor tends to smooth out the charge density. 

The results of HF theory and Hartree theory are very similar because the exchange interaction between orbitals on 
different sites is small. The orbitals obtained within unrestricted HF theory (and Hartree theory) are very strongly 
localized and the kinetic energy within HF theory is larger than in our DMC calculations with the fully optimized 
trial wave functions. We have defined the correlation energy to be the difference between the exact and unrestricted 
HF energies. With this definition, and in the low density limit, Wigner crystals are weakly correlated systems. The 
inclusion of correlation in a Wigner crystal wave function beyond the unrestricted HF level results in the electronic 
charge density spreading out from the lattice sites. In this sense correlation delocalizes the electrons. Although the 
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FIG. 8: Spin-independent PCF for the antiferromagnetic floating crystal and unpolarized fluid at r s — 110. Solid line: 
extrapolated estimator for the crystal; dotted line: HF result; dashed-dotted line: VMC result obtained using the DMC 
optimized orbitals but without a Jastrow factor. The dashed line shows the extrapolated estimator for the unpolarized fluid. 
The statistical errors are less than 0.003 except for the dashed-dotted line around r = 0. 



correlation energies of the crystal phases are smaller than those of the fluid phases at the same density, the use of HF 
orbitals within the trial wave functions results in larger fixed node errors for the crystal phases. 

The variational energy for a floating Wigner crystal is lower than that of the fixed crystal by the kinetic energy of 
the center of mass, which is a negligible energy per particle for large systems. The expectation value of any operator 
that depends only on the difference between electron coordinates is the same for the fixed and floating crystals. We 
can therefore obtain the PCFs of a floating Wigner crystal rather simply from calculations on the fixed crystal. 
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